* Gary LaMotte's Mathematica solutions to exercises from the 2012 Gravitational Wave Astronomy Summer School

*DAY 3*

*1.  &     2.  ----  Done on DAY 2*

*3.  MATCHED FILTERING IN THE FREQUENCY DOMAIN*

*a.*

samprate=16384;

data =  Flatten[Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day2data.txt", "Table"]];

datasecs = Length[data]/samprate

tdata =  Table[(r-1.0)/samprate  , {r,Length[data]}];

fftdata = Fourier[data];

*b.*

template= Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day2template.txt", "List"];

exttemplate = PadRight[template, 262144];

ffttemplate = Fourier[exttemplate ];

*c.*

multiplybins = fftdata*(ffttemplate)\[Conjugate];

freqcorr = InverseFourier[multiplybins];

*d.*

filtoutputtimeseries = Transpose[{tdata, freqcorr}];

ListLinePlot[filtoutputtimeseries , PlotLabel -> NonOptimal Output ,PlotRange-> {{0,16}, {-.00000000000000000012,.00000000000000000021}}, AxesOrigin->{-.00000000000000006,0},  Frame->True, AxesLabel ->  {Time (Secs), Amplitude}]

*The graph has the same shape as the time domain matched filtering from Day 2, but has different actual values (differing only by a proportionality constant).  
Frequency domain matched filtering is much faster to compute than time domain matched filtering.*

*4.  OPTIMAL MATCHED FILTERING USING FREQUENCY WEIGHTING*

*a.*

noise =  Flatten[Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day3asd.txt", "Table"]];

noisespectrum = noise[[1;;262144]];

noisepsd = Abs[(noisespectrum )^2];

*b.*

multiplybinsdividebypsd = fftdata*(ffttemplate)\[Conjugate]/noisepsd;

optcorr = InverseFourier[multiplybinsdividebypsd ];

*c.*

optimalmatchtimeseries = Transpose[{tdata, optcorr}];

ListLinePlot[ optimalmatchtimeseries, PlotLabel -> Optimal Output,PlotRange-> {{0,16}, {-.0000000000000000027,.000000000000000006}}, AxesOrigin->{-.00000000000000006,0},  Frame->True, AxesLabel ->  {Time (Secs), Amplitude}]

*Optimal correlation is only slightly different from non-optimal correlation.  This is because the noise power spectrum does not vary much in the limited frequency range of the template from appox 35 to 250Hz.  SNR is approx 14%  improved.*

*d.*

firstsecsnonoptimal =  freqcorr  [[1;;50000]];

snrnonoptimal = Max[freqcorr ]/StandardDeviation[firstsecsnonoptimal]

firstsecsoptcorr =  optcorr [[1;;50000]];

snroptimal = Max[optcorr ]/StandardDeviation[firstsecsoptcorr ]

snrimprovement = snroptimal/snrnonoptimal

*e.*

periodo = Abs[(fftdata)^2];

periodoproduct=fftdata*(ffttemplate)\[Conjugate]/periodo;

pcorr = InverseFourier[periodoproduct];

firstpcorr = pcorr[[1;;50000]];

periodosnr =Max[pcorr]/StandardDeviation[firstpcorr]

*Next to smooth the periodo a convolution with a moving average filter (ker)will be performed.*

ker = PadRight[{1,1},25,1]

speriodo =ListConvolve[ker/Length[ker],periodo];

reciprocalsperiodo = 1/speriodo;

sproduct=fftdata*(ffttemplate)\[Conjugate]*PadRight[reciprocalsperiodo,262144];

spcorr = Re[InverseFourier[sproduct]];

firstspcorr = spcorr[[1;;50000]];

speriodosnr =Max[spcorr]/StandardDeviation[firstspcorr]

*Smoothing the PSD gives a higher SNR vs non-smoothed.  Smoothed PSD still does not have as high a SNR as does optimal matched filtering (84 vs 92).*

*Next a plot of the analytic PSD and the smoothed periodogram overlaid.*

psdfreq = Table[(p-1.0)/datasecs, {p,Length[noisepsd]}];

PSDperiodogram = Transpose[{psdfreq, noisepsd  }];

speriodoplot=PadRight[speriodo, 262144,speriodo[[Length[speriodo]]]];

speriodoperiodogram = Transpose[{psdfreq, speriodoplot  }];

ListLogPlot[{speriodoperiodogram,PSDperiodogram},PlotRange->{{0,512},{.01, 1000}},Joined->True, PlotLabel->ComparisonPeriodogram, AxesLabel-> {Frequency, Amp^2}]

*Smoothed Periodogram vs. analytic PSD:  
The signal (blue)demonstrates increased power from 40 to 230 Hz above the noise level (red).  This is due to the signal in the data.*